/**
 * @file ebex_de.h
 *
 * @date Jan 12, 2010
 * @author seth
 *
 * @brief This file is part of EBEX State Estimator, created for the EBEX project
 *
 * This software is copyright (C) 2010 Columbia University
 *
 * EBEX State Estimator is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * EBEX State Estimator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with EBEX State Estimator; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * @note Differential evolution was conceived by Ken Price and Rainer Storn.
 * For details see http://www.icsi.berkeley.edu/~storn/code.html
 */

#include <ebex_genetic.h>
#include <ebex_quaternion.h>
#include <ebex_sort.h>
#include <ebex_genetic_stats.h>
#include <ebex_de.h>

/**
 * Sets the differential evolution parameters for the population
 * @param m_pop Population pointer
 * @param m_strategy Differential Evolution strategy (#diff_evo_strategy_t)
 * @param m_crossover_type Type of crossover (binomial/exponential)
 * @param m_numvec Number of vectors to use for the differencing
 * @param m_minweight Minimum weight to use for crossover chromosomes
 * @param m_maxweight Maximumum weight to use for crossover chromosomes
 * @param m_crossover_ratio
 */
void ebex_diffevolv_setparameters(population_t *m_pop, const diff_evo_strategy_t m_strategy,
		const diff_evo_crossover_t m_crossover_type, const int m_numvec, const double m_minweight,
		const double m_maxweight, const double m_crossover_ratio, const double m_opposition_rate,
		const diff_evo_reference_t m_reference, const bool m_scatter)
{

	if (!m_pop)
		die("Null pointer to population structure passed.");

	elog( log_level_verbose, "Population's differential evolution parameters set");

	if (m_pop->de_params == NULL)
	{
		if (!(m_pop->de_params = e_malloc(sizeof(diff_evo_t))))
			die("Unable to allocate memory");
	}

	m_pop->de_params->strategy = m_strategy;
	m_pop->de_params->crossover_type = m_crossover_type;
	m_pop->de_params->num_dv = m_numvec;
	m_pop->de_params->min_weight = m_minweight;
	m_pop->de_params->max_weight = m_maxweight;
	m_pop->de_params->crossover_ratio = m_crossover_ratio;
	m_pop->de_params->opposition_rate = m_opposition_rate;
	m_pop->de_params->reference = m_reference;
	m_pop->de_params->scatter = m_scatter;
	m_pop->de_params->scatter_chromosome = 0;

	if (m_pop->de_params->crossover_ratio < 0.0 || m_pop->de_params->crossover_ratio > 1.0)
		dief("Invalid Differential Evolution Crossover Ratio %f.  Must be (0,1]", m_pop->de_params->crossover_ratio);

	if (m_pop->de_params->min_weight > m_pop->de_params->max_weight)
	{
		elog(log_level_warning, "Differential Evolution minimum weight less than max weight.  Inverting values");
		EBEX_SWAP(m_pop->de_params->min_weight, m_pop->de_params->max_weight);
	}
	return;
}


static inline void ebex_diffevolv_select_rnd(int *m_permutation, int m_cnt, int m_size, int m_skip)
{
	int pos;

	for (int i = 0; i < m_cnt; i++)
	{
		do
		{
			pos = i + ran_uint_uniform(m_size - i);
		} while (m_permutation[pos] == m_skip);

		EBEX_SWAP(m_permutation[i], m_permutation[pos]);
	}

	return;
}

static inline void ebex_diffevolv_crossover_seed_double(population_t *m_pop, organism_t *m_entity, int *m_permutation,
                                                        double m_weight, int m_seed, size_t m_chromosome, size_t m_allele)
{
	double diff = 0.0;
	size_t num = m_pop->de_params->num_dv;

	for (size_t i = 0; i < num; i++)
	{
		diff += (m_pop->sorted_organisms[m_permutation[i]]->chromosome[m_chromosome].dbl_allele[m_allele] -
				m_pop->sorted_organisms[m_permutation[(num * 2 - 1) - i]]->chromosome[m_chromosome].dbl_allele[m_allele]);
	}

	m_entity->chromosome[m_chromosome].dbl_allele[m_allele] = m_pop->sorted_organisms[m_seed]->chromosome[m_chromosome].dbl_allele[m_allele] + m_weight * diff;
}

static inline void ebex_diffevolv_crossover_seed_quaternion(population_t *m_pop, organism_t *m_entity, int *m_permutation,
                                                            double m_weight, int m_seed, size_t m_chromosome, size_t m_allele)
{
	size_t num = m_pop->de_params->num_dv;
	ebex_quat_t diff;
	ebex_quat_t total_diff;

	ebex_quat_error(&total_diff, &m_pop->sorted_organisms[m_permutation[0]]->chromosome[m_chromosome].q_allele[m_allele],
					&m_pop->sorted_organisms[m_permutation[(num * 2 - 1)]]->chromosome[m_chromosome].q_allele[m_allele]);

	for (size_t i = 1; i < num; i++)
	{
		ebex_quat_error(&diff, &m_pop->sorted_organisms[m_permutation[i]]->chromosome[m_chromosome].q_allele[m_allele],
				&m_pop->sorted_organisms[m_permutation[(num * 2 - 1) - i]]->chromosome[m_chromosome].q_allele[m_allele]);
		ebex_quat_mul(&total_diff, &total_diff, &diff);
	}

	ebex_quat_scale_angle(&total_diff, &total_diff, m_weight);
	ebex_quat_mul(&m_entity->chromosome[m_chromosome].q_allele[m_allele], &m_pop->sorted_organisms[m_seed]->chromosome[m_chromosome].q_allele[m_allele], &total_diff);
}

static inline void ebex_diffevolv_crossover_randbest_double(population_t *m_pop, organism_t *m_entity, int *m_permutation,
                                                            double m_weight, int m_seed, size_t m_chromosome, size_t m_allele)
{
	double diff;
	size_t num = m_pop->de_params->num_dv;

	diff = m_pop->sorted_organisms[m_seed]->chromosome[m_chromosome].dbl_allele[m_allele] - m_entity->chromosome[m_chromosome].dbl_allele[m_allele];
	for (size_t i = 0; i < num; i++)
	{
		diff += (m_pop->sorted_organisms[m_permutation[i]]->chromosome[m_chromosome].dbl_allele[m_allele] -
				m_pop->sorted_organisms[m_permutation[(num * 2 - 1) - i]]->chromosome[m_chromosome].dbl_allele[m_allele]);
	}

	m_entity->chromosome[m_chromosome].dbl_allele[m_allele] += m_weight * diff;
}

static inline void ebex_diffevolv_crossover_randbest_quaternion(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight, int m_seed, size_t m_chromosome, size_t m_allele)
{
	size_t num = m_pop->de_params->num_dv;
	ebex_quat_t diff;
	ebex_quat_t total_diff = {0.0};

	ebex_quat_error(&total_diff, &m_pop->sorted_organisms[m_seed]->chromosome[m_chromosome].q_allele[m_allele],
					&m_entity->chromosome[m_chromosome].q_allele[m_allele]);

	for (size_t i = 0; i < num; i++)
	{
		ebex_quat_error(&diff, &m_pop->sorted_organisms[m_permutation[i]]->chromosome[m_chromosome].q_allele[m_allele],
				&m_pop->sorted_organisms[m_permutation[(num * 2 - 1) - i]]->chromosome[m_chromosome].q_allele[m_allele]);
		ebex_quat_mul(&total_diff, &total_diff, &diff);
	}

	ebex_quat_scale_angle(&total_diff, &total_diff, m_weight);
	ebex_quat_mul(&m_entity->chromosome[m_chromosome].q_allele[m_allele], &m_pop->sorted_organisms[m_permutation[0]]->chromosome[m_chromosome].q_allele[m_allele], &total_diff);
}

static void ebex_diffevolv_crossover_bin_best(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight, int m_best)
{
	size_t allele;
	size_t chromosome;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (chromosome = 0; chromosome < m_pop->num_chromosomes; chromosome++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_seed_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_seed_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
		{
			if (ran_bool_prob(m_pop->de_params->crossover_ratio)) crossover(m_pop, m_entity, m_permutation, m_weight, m_best, chromosome, allele);
		}
	}

}

static void ebex_diffevolv_crossover_bin_rand(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight)
{
	size_t initial_chromosome = ran_uint_uniform(m_pop->num_chromosomes);
	size_t initial_allele = ran_uint_uniform(m_entity->chromosome[initial_chromosome].length);
	size_t allele = initial_allele;
	size_t chromosome = initial_chromosome;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (chromosome = 0; chromosome < m_pop->num_chromosomes; chromosome++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_seed_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_seed_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
		{
			if (ran_bool_prob(m_pop->de_params->crossover_ratio)) crossover(m_pop, m_entity, m_permutation, m_weight, m_permutation[m_pop->de_params->num_dv], chromosome, allele);
		}
	}
}

static void ebex_diffevolv_crossover_bin_randtobest(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight, int m_best)
{
	size_t allele;
	size_t chromosome;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (chromosome = 0; chromosome < m_pop->num_chromosomes; chromosome++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_randbest_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_randbest_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
		{
			if (ran_bool_prob(m_pop->de_params->crossover_ratio)) crossover(m_pop, m_entity, m_permutation, m_weight, m_permutation[m_pop->de_params->num_dv], chromosome, allele);

		}
	}
}

static void ebex_diffevolv_crossover_exp_best(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight, int m_best)
{
	size_t initial_chromosome = ran_uint_uniform(m_pop->num_chromosomes);
	size_t initial_allele = ran_uint_uniform(m_entity->chromosome[initial_chromosome].length);
	size_t allele = initial_allele;
	size_t chromosome = initial_chromosome;
	size_t i, j;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (i = 0; i < m_pop->num_chromosomes; i++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_seed_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_seed_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		allele %= m_entity->chromosome[chromosome].length;
		for (j = 0; j < m_entity->chromosome[chromosome].length; j++)
		{
			if (!ran_bool_prob(m_pop->de_params->crossover_ratio)) return;

			crossover(m_pop, m_entity, m_permutation, m_weight, m_best, chromosome, allele);
			allele = (allele + 1) % m_entity->chromosome[chromosome].length;

		}
		chromosome = (chromosome + 1) % m_pop->num_chromosomes;
	}
}

static void ebex_diffevolv_crossover_exp_rand(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight)
{
	size_t initial_chromosome = ran_uint_uniform(m_pop->num_chromosomes);
	size_t initial_allele = ran_uint_uniform(m_entity->chromosome[initial_chromosome].length);
	size_t allele = initial_allele;
	size_t chromosome = initial_chromosome;
	size_t i,j;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (i = 0; i < m_pop->num_chromosomes; i++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_seed_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_seed_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		allele %= m_entity->chromosome[chromosome].length;
		for (j = 0; j < m_entity->chromosome[chromosome].length; j++)
		{
			if (!ran_bool_prob(m_pop->de_params->crossover_ratio)) return;

			crossover(m_pop, m_entity, m_permutation, m_weight, m_permutation[m_pop->de_params->num_dv], chromosome, allele);
			allele = (allele + 1) % m_entity->chromosome[chromosome].length;
		}
		chromosome = (chromosome + 1) % m_pop->num_chromosomes;
	}
}

static void ebex_diffevolv_crossover_exp_randtobest(population_t *m_pop, organism_t *m_entity, int *m_permutation, double m_weight, int m_best)
{
	size_t initial_chromosome = ran_uint_uniform(m_pop->num_chromosomes);
	size_t initial_allele = ran_uint_uniform(m_entity->chromosome[initial_chromosome].length);
	size_t allele = initial_allele;
	size_t chromosome = initial_chromosome;
	size_t i,j;

	void (*crossover)(population_t*,organism_t*,int*,double,int,size_t,size_t);

	for (i = 0; i < m_pop->num_chromosomes; i++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				crossover = ebex_diffevolv_crossover_randbest_double;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				crossover = ebex_diffevolv_crossover_randbest_quaternion;
				break;
			default:
				die("Invalid crossover type for Differential Evolution");
				break;
		}

		allele %= m_entity->chromosome[chromosome].length;
		for (j = 0; j < m_entity->chromosome[chromosome].length; j++)
		{
			if (!ran_bool_prob(m_pop->de_params->crossover_ratio)) return;

			crossover(m_pop, m_entity, m_permutation, m_weight, m_best, chromosome, allele);
			allele = (allele + 1) % m_entity->chromosome[chromosome].length;
		}
		chromosome = (chromosome + 1) % m_pop->num_chromosomes;
	}
}

static inline void ebex_diffevolv_double_opposite(organism_t *m_entity, size_t m_chromosome, size_t m_allele, const organism_t *m_reference)
{
	double diff = m_reference->chromosome[m_chromosome].dbl_allele[m_allele] - m_entity->chromosome[m_chromosome].dbl_allele[m_allele];

	m_entity->chromosome[m_chromosome].dbl_allele[m_allele] += (diff + diff);
}

static inline void ebex_diffevolv_quaternion_opposite(organism_t *m_entity, size_t m_chromosome, size_t m_allele, const organism_t *m_reference)
{
	ebex_quat_t diff;

	ebex_quat_error(&diff, &(m_entity->chromosome[m_chromosome].q_allele[m_allele]),
	                &(m_reference->chromosome[m_chromosome].q_allele[m_allele]));
	ebex_quat_scale_angle(&diff, &diff, 2.0);

	ebex_quat_mul(&(m_entity->chromosome[m_chromosome].q_allele[m_allele]), &(m_entity->chromosome[m_chromosome].q_allele[m_allele]), &diff);

}

static void ebex_diffevolv_generate_opposite(population_t *m_pop, organism_t *m_entity, organism_t *m_reference)
{
	size_t allele;
	size_t chromosome;

	void (*opposite)(organism_t*,size_t,size_t,const organism_t*);

	for (chromosome = 0; chromosome < m_pop->num_chromosomes; chromosome++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				opposite = ebex_diffevolv_double_opposite;
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				opposite = ebex_diffevolv_quaternion_opposite;
				break;
			default:
				die("Invalid chromosome type for Differential Evolution");
				break;
		}

		for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
		{
			opposite(m_entity, chromosome, allele, m_reference);
		}
	}
}

static void ebex_diffevolv_generate_mean(population_t *m_pop, organism_t *m_entity)
{
	size_t allele;
	size_t chromosome;

	for (chromosome = 0; chromosome < m_pop->num_chromosomes; chromosome++)
	{
		switch(m_entity->chromosome[chromosome].type)
		{
			case CHROMOSOME_TYPE_DOUBLE:
				for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
				{
					m_entity->chromosome[chromosome].dbl_allele[allele] = genetic_allele_mean_double(m_pop, chromosome, allele, m_pop->orig_size);
				}
				break;
			case CHROMOSOME_TYPE_QUATERNION:
				for (allele = 0; allele < m_entity->chromosome[chromosome].length; allele++)
				{
					genetic_allele_mean_quaternion(m_pop, chromosome, allele, &(m_entity->chromosome[chromosome].q_allele[allele]), m_pop->orig_size);
				}
				break;
			default:
				die("Invalid chromosome type for Differential Evolution");
				break;
		}
	}
}

static void ebex_differential_evolution_opposition(population_t *pop)
{
	size_t i;
	organism_t *new_entity;

	pop->orig_size = pop->size;
	organism_t *reference_organism;

	switch(pop->de_params->reference)
	{
		case DIFF_EVO_REFERENCE_MEAN:
			reference_organism = genetic_get_new_organism(pop);
			ebex_diffevolv_generate_mean(pop, reference_organism);
			pop->evaluate(pop, reference_organism);
			break;
		case DIFF_EVO_REFERENCE_BEST:
			reference_organism = pop->sorted_organisms[0];
			break;
		case DIFF_EVO_REFERENCE_RANDOM:
			reference_organism = pop->sorted_organisms[ran_uint_uniform(pop->size)];
			break;
		default:
			die("Invalid differential evolution reference type!");
			break;
	}

#	pragma omp parallel for 	\
		shared(pop, reference_organism) private(i) schedule(static)
	for (i = 0; i < pop->orig_size; i++)
	{
		new_entity = genetic_organism_clone(pop, pop->sorted_organisms[i]);
		ebex_diffevolv_generate_opposite(pop, new_entity, reference_organism);

		if (pop->verify(pop, new_entity))
		{
			pop->evaluate(pop, new_entity);
		}
	}

	sort_population(pop);
	genetic_genocide(pop, pop->stable_size);
}

void ebex_diffevolv_evolve_organism(population_t *m_pop, organism_t *m_entity, int *m_permutation_array, double m_weight, int m_ref, int m_skip)
{
	if (m_pop->de_params->strategy == DIFF_EVO_STRATEGY_RANDOM)
	{
		ebex_diffevolv_select_rnd(m_permutation_array, 2 * m_pop->de_params->num_dv + 1, m_pop->orig_size, m_skip);
	}
	else
	{
		ebex_diffevolv_select_rnd(m_permutation_array, 2 * m_pop->de_params->num_dv, m_pop->orig_size, m_skip);
	}

	switch (m_pop->de_params->crossover_type)
	{
		case DIFF_EVO_CROSSOVER_BINOMIAL:

			switch(m_pop->de_params->strategy)
			{
				case DIFF_EVO_STRATEGY_BEST:
					ebex_diffevolv_crossover_bin_best(m_pop, m_entity, m_permutation_array, m_weight, m_ref);
					break;
				case DIFF_EVO_STRATEGY_RANDOM:
					ebex_diffevolv_crossover_bin_rand(m_pop, m_entity, m_permutation_array, m_weight);
					break;
				case DIFF_EVO_STRATEGY_RANDOMTOBEST:
					ebex_diffevolv_crossover_bin_randtobest(m_pop, m_entity, m_permutation_array, m_weight, m_ref);
					break;
				default:
					die("Unknown differential evolution strategy.");
					break;
			}
			break;
		case DIFF_EVO_CROSSOVER_EXPONENTIAL:
			switch(m_pop->de_params->strategy)
			{
				case DIFF_EVO_STRATEGY_BEST:
					ebex_diffevolv_crossover_exp_best(m_pop, m_entity, m_permutation_array, m_weight, m_ref);
					break;
				case DIFF_EVO_STRATEGY_RANDOM:
					ebex_diffevolv_crossover_exp_rand(m_pop, m_entity, m_permutation_array, m_weight);
					break;
				case DIFF_EVO_STRATEGY_RANDOMTOBEST:
					ebex_diffevolv_crossover_exp_randtobest(m_pop, m_entity, m_permutation_array, m_weight, m_ref);
					break;
				default:
					die("Unknown differential evolution strategy.");
					break;
			}
			break;
		default:
			die("Unknown differential evolution crossover method");
			break;
	}
}

size_t ebex_differential_evolution(population_t *m_pop, const size_t max_generations)
{
	size_t generation = 0;
	size_t i;
	size_t best;
	int *permutation;
	organism_t *new_entity;
	double weight;

	if (!m_pop)
		die("NULL pointer to population structure passed.");
	if (!m_pop->de_params)
		die("Missing parameters to differential evolution");
	if (!m_pop->evaluate)
		die("Population's evaluation callback is undefined.");
	if (!m_pop->verify)
		die("Population's verification callback is undefined.");
	if (!m_pop->chromosome_constructor)
		die("Population's chromosome constructor callback is undefined.");
	if (!m_pop->chromosome_destructor)
		die("Population's chromosome destructor callback is undefined.");
	if (m_pop->stable_size < 6)
		die("Population's stable size is too small.  (Must be at least 6)");
	if (m_pop->de_params->crossover_ratio < 0.0 || m_pop->de_params->crossover_ratio > 1.0)
		die("Invalid crossover_ratio.");

	elog(log_level_verbose, "The differential evolution has begun!");

	m_pop->generation = 0;

	/*
	 * Score the initial population members.
	 */
	if (m_pop->size < m_pop->stable_size)
		genetic_fill_pop(m_pop, m_pop->stable_size - m_pop->size);

	if (m_pop->sorted_organisms[0]->fitness <= GENETIC_MIN_FITNESS)
		m_pop->evaluate(m_pop, m_pop->sorted_organisms[0]);

#	pragma omp parallel for \
	   shared(m_pop) private(i) \
	   schedule(static)
	for (i = 0; i < m_pop->size; i++)
	{
		if (m_pop->sorted_organisms[i]->fitness <= GENETIC_MIN_FITNESS)
		{
			m_pop->evaluate(m_pop, m_pop->sorted_organisms[i]);
		}
	}

	sort_population(m_pop);

	/*
	 * Prepare arrays to store permutations.
	 */
	if (!(permutation = e_malloc(sizeof(int)*m_pop->size)))
		die("Unable to allocate memory");

	for (i = 0; i < m_pop->size; i++)
		permutation[i] = i;

	/*
	 * Do all the generations:
	 *
	 * Stop when (a) max_generations reached, or
	 *           (b) "pop->generation_hook" returns false.
	 */
	while (generation < max_generations)
	{
		m_pop->generation = ++generation;
		m_pop->orig_size = m_pop->size;

		if (m_pop->pre_generation_hook && (!m_pop->pre_generation_hook(generation, m_pop, NULL))) break;

		/// Dither the weighting vector if requested (min_weight < max_weight
		weight = ran_double_range_uniform(m_pop->de_params->min_weight, m_pop->de_params->max_weight);
		
		/// Find the best solution.
		best = 0;

#		pragma omp parallel for 	\
		   shared(m_pop) private(i) schedule(static)
		for (i = 0; i < m_pop->orig_size; i++)
		{

			new_entity = genetic_organism_clone(m_pop, m_pop->sorted_organisms[i]);

			ebex_diffevolv_evolve_organism(m_pop, new_entity, permutation, weight, best, i);

			/**
			 * Evaluate our newest solution.  If it is valid, evaluate and give it a fitness value,
			 * otherwise, set its fitness to the minimum, so that it will be culled during the genocide
			 */
			if (m_pop->verify(m_pop, new_entity))
			{
				m_pop->evaluate(m_pop, new_entity);
			}
			else
			{
				new_entity->fitness = GENETIC_MIN_FITNESS;
			}
		}

		/**
		 * After creating our new individuals with new fitnesses, we sort and kill off the least fit
		 * members of our population.  These may be parents or offspring.
		 */
		sort_population(m_pop);
		genetic_genocide(m_pop, m_pop->stable_size);


		if (m_pop->de_params->reference != DIFF_EVO_REFERENCE_NONE && ran_unit_uniform() < m_pop->de_params->opposition_rate)
		{
			ebex_differential_evolution_opposition(m_pop);
		}

		if (m_pop->post_generation_hook && (!m_pop->post_generation_hook(generation, m_pop, NULL))) break;

	} /// Generation loop

	e_free(permutation);

	return generation;
}

